#include "JaccardIndex.h"

JaccardIndex::JaccardIndex()
{
}

JaccardIndex::~JaccardIndex()
{
}

double** JaccardIndex::calc(map<string,double*>& Abundance, int SampleNum, double** JaccardIndexMatrix)
{
    int tPairs=SampleNum*(SampleNum-1)/2;
    int* tX=new int[tPairs];         //sum of OTU of first community
    memset(tX,0,sizeof(int)*tPairs);  
    int* tY=new int[tPairs];         //sum of OTU of second community
    memset(tY,0,sizeof(int)*tPairs);
    int* tC=new int[tPairs];
    memset(tC,0,sizeof(int)*tPairs);
    map<string,double*>::iterator it=Abundance.begin();
    while(it!=Abundance.end())
    {
        for(int i=0;i<SampleNum-1;i++)
            for(int j=i+1;j<SampleNum;j++)
            {
                int tIndex=i*SampleNum-i*(i+1)/2+j-i-1;
                if((*it).second[i]!=0)
                    tX[tIndex]++;
                if((*it).second[j]!=0)
                    tY[tIndex]++;
                if((*it).second[i]!=0&&(*it).second[j]!=0)
                    tC[tIndex]++;
            }
        it++;
    }
    for(int i=0;i<SampleNum;i++)
        for(int j=0;j<SampleNum;j++)
        {
            if(i==j)
            {
                JaccardIndexMatrix[i][j]=0;
            }
            if(i>j)  
                JaccardIndexMatrix[i][j]=JaccardIndexMatrix[j][i];
            if(i<j)
            {
                int tIndex=i*SampleNum-i*(i+1)/2+j-i-1;
                JaccardIndexMatrix[i][j]=(double)tC[tIndex]/(tX[tIndex]+tY[tIndex]-tC[tIndex]);
            }
            //cout<<JaccardIndexMatrix[i][j];
            //if(j==SampleNum-1)
                //cout<<"\n";
            //else
                //cout<<"\t";
        }
    delete[] tX;
    delete[] tY;
    delete[] tC;
    return JaccardIndexMatrix;
}
